{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### <span style=\"color:#cc7a00\">Critical Zone Observatory - Modeling Institute 2016</span>\n",
    "\n",
    "<img style=\"float: right;\" src=\"http://hydrocomplexity.net/images/dharaD.png\" width=\"70px\" hspace=5 /> <br/>\n",
    "\n",
    "# Tutorial: Introduction to Dhara<sup>$\\dagger$</sup> model\n",
    "*<sup>$\\dagger$</sup>Sanskrit, noun: earth, flow <br>*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. What is Dhara?\n",
    "\n",
    "[Dhara](https://hydrocomplexity.github.io/Dhara/) is an open-source, high-performance computing framework for modeling complex processes in the __Critical Zone__ (CZ). The model is designed to be modular in structure with the aim to create uniform, efficient, and interdisciplinary tools to facilitate and leverage process modeling. It also provides flexibility to maintain, collaborate, and co-develop new components by the scientific community. Currently, Dhara consists of two main components (See Figure 1):\n",
    "* A <span style=\"color:#cc7a00\">**multi-layer canopy**</span> model to simulate within-canopy processes in the aboveground,\n",
    "* An <span style=\"color:#cc7a00\">**integrated flow**</span> model which includes overland and subsurface modules to simulate vertical and horizontal fluxes on the surface and in the belowground system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://hydrocomplexity.net/images/dhara/Dhara_main.png\" alt=\"\" width=\"600\">\n",
    "\n",
    "</br>\n",
    "<center>*__Figure 1__. Schematic of Dhara model on hybrid CPU-GPU parallel computing architecture.*</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dhara is primarily developed to utilize *high-resolution topographic* data (particularly *lidar*) for simulations of complex (i.e. ecohydrologic and biogeochemical) processes dominated by topographic *controls* and *variability* across scales in the CZ. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Organization\n",
    "\n",
    "Dhara modeling framework uses hybrid parallel computing which combines the advantages of [Message Passing Interface (MPI)](https://www.open-mpi.org/) and [GPU parallel computing](http://www.nvidia.com/object/what-is-gpu-computing.html) for optimal performance on large-scale simulations. \n",
    "The model is written using [CUDA C/C++](https://blogs.nvidia.com/blog/2012/09/10/what-is-cuda-2/) parallel computing platform and programming model that are based on the industry-standard C/C++.\n",
    "With this approach, users can accelerate their code by moving the computationally intensive, data parallelism portion to a GPU and still utilize MPI for portions that are best suited to task parallelism for optimal performances."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Structure\n",
    "The structure of Dhara source code is arranged in 4 different levels as shown in Figure 2.\n",
    "* __Level 1__ - main program entry point (`main.cc`): Program starts and controls MPI environments and numerical modules;\n",
    "* __Level 2__ - numerical (`numerical.cc`) and MPI (`mpienv.cc`) modules: While `mpienv.cc` provides functions for MPI , `numerical.cc` includes functions and links to level 3 files for inputs/outputs (I/O) and numerical simulations.\n",
    "* __Level 3__ - processing modules: includes files that provide tools to process I/O (`io.cc`) in `Dhara`, parse configurations (`parser.cc`), and execute host (`host.cc`) and device (`device.cu`) functions.\n",
    "* __Level 4__ - model solvers: includes source files for numerical solvers for different modules used in Dhara. Most of the level 4 files are linked with the `device.cu` where all the models communicate to each other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://hydrocomplexity.net/images/dhara/dhara_codestructure.png?nocache=<?php echo time(); ?>\" alt=\"\" width=\"600\">\n",
    "\n",
    "</br>\n",
    "<center>*__Figure 2__. Organization of source files in Dhara. User-defined classes for Dhara model is shown on the right*</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "User-defined classes are defined in header files. Dhara source code is built and compiled using `Makefile`:\n",
    "* All `*.cc` files are compiled with `GNU gcc` compiler with a link to CUDA (`lcudart`) library.\n",
    "* All `*.cu` files are compiled with `nvcc` compiler with a link to MPI. \n",
    "* NetCDF library is linked to Dhara for I/O using array-oriented scientific data\n",
    "* Eigen C++ template library is used for linear algebra in MLCan model\n",
    "* CUSP library is used for solving sparse linear systems in GPU for integrated flow modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Workflow\n",
    "Figure 3 presents the workflow of Dhara model in a GPU computing node on supercomputers. \n",
    "Typically, GPU computing nodes have a GPU device and a number of CPUs.\n",
    "When computing in parallel, a number of processes run the model at the same time. The first MPI process, which we call the *master* or *root* process, also controls data movement between processes and communication with the GPU device."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://hydrocomplexity.net/images/dhara/dhara_blocks.png?cache=none\" alt=\"\" width=\"700\">\n",
    "\n",
    "</br>\n",
    "<left>*__Figure 3. Workflow of Dhara model in hybrid CPU-GPU parallel computing architecture__ - MPI is activated at the beginning of the program for parallel computing in MLCan model. Only the master process (rank=0) initializes CUDA environment for integrated flow modeling. Other processes can see the device (GPU) but don't need to work on it. Simulations of MLCan model on all processes are gathered/scattered to/from the master process for communicating with device for integrated flow modeling. After simulations are completed, MPI and CUDA environments are terminated.*</left>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown in the figure, a number of input data and parameters is required to run the model. The list of these inputs for `Dhara` are shown below:\n",
    "\n",
    "1. **Multi-layer canopy model**\n",
    "    * *Vegetation* (each plant species)\n",
    "        * Canopy: Leaf Area Index - `LAI`, Leaf Area Density - `LAD`\n",
    "        * Root: root density - `rootfr`\n",
    "        * Parameters: Plant specie parameters are also required for running MLCan module. A list of important parameters for MLCan module can be found in its original document [here](http://onlinelibrary.wiley.com/doi/10.1029/2010JG001340/full).\n",
    "        \n",
    "    * *Forcings*\n",
    "        * Precipiation\n",
    "        * Air temperature\n",
    "        * Incoming shortwave radiation\n",
    "        * Longwave radiation (optional)\n",
    "        * Vapor pressure\n",
    "        * Atmospheric pressure\n",
    "        * Wind velocity\n",
    "        * Sun zenith angle\n",
    "        \n",
    "2. **Integrated flow model**\n",
    "    * *Overland flow*\n",
    "        * Topographic data\n",
    "        * Initial and boundary conditions\n",
    "        * Parameters: Manning's coefficients\n",
    "    * *Subsurface flow*\n",
    "        * Soil parameters\n",
    "        * Initial and boundary conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Getting started\n",
    "\n",
    "### 3.1 Bioenergy crop example\n",
    "The study site in this example is located at the energy farm managed by [Energy Biosciences Institute (EBI)](http://www.energybiosciencesinstitute.org) at the [University of Illinois](https://illinois.edu). \n",
    "We will refer to the site as EBI from now on. \n",
    "Figure 4 shows a satellite image of EBI site taken in 2008."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://hydrocomplexity.net/images/dhara/ebi_map.png?cache=none\" alt=\"\" width=\"800\">\n",
    "\n",
    "</br>\n",
    "<center>*__Figure 4. Satellite image of EBI study site__ - The sites is 400m x 400m and has 4 subplots that plant several energy crops including corn (2 plots, top), miscanthus (1 plot, bottom left), and switchgrass (1 plot, bottom right).*</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to make sure that all inputs are valid and reasonable before running Dhara model. \n",
    "First, let's take a look and visualize all input data at this site using jupyter notebook. \n",
    "Inputs are located at: \n",
    "\n",
    "<pre>\n",
    "czo-mi$ Dhara_folder/examples/1-EBI\n",
    "</pre>    \n",
    "\n",
    "Similar to the batch environment in Unix, we can set Dhara folder as a variable named `$Dhara` in this jupyter notebook for fast changing directory in the future.\n",
    "\n",
    "<pre>\n",
    "Dhara='~/czomi/users/<span style=\"color:#cc7a00\"><b>YOUR_USERNAME</b></span>/Dhara'\n",
    "</pre>\n",
    "\n",
    "Change directory to Dhara_folder:\n",
    "<pre>\n",
    "cd $Dhara\n",
    "</pre>\n",
    "\n",
    "Now, list the content in Dhara folder:\n",
    "<pre>\n",
    "ls -l\n",
    "</pre>\n",
    "\n",
    "Then, list the content in the EBI site example to screen:\n",
    "<pre>\n",
    "ls -l examples/1-EBI/\n",
    "</pre>\n",
    "\n",
    "\n",
    "<sub> __Note__: these above commands must be executed in separate jupyter python cells to be effective.</sub>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__TODO:__ create new cells and execute above commands"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see several __configuration__ (`.cfg`) and __netcdf__ (`.nc`) files in the folder. \n",
    "They all contain information needed to drive the model.\n",
    "While configuration files includes *constant parameters* and *model options*, netcdf files often includes *temporal* and *spatial* datasets$^{\\ddagger}$.\n",
    "There is one <span style=\"color:#cc7a00\"> __portable batch system__</span> (`.pbs`) file for submitting jobs on supercomputer that will be explained later.\n",
    "We'll use [Python](https://www.python.org) to visualize and analyze all the inputs and model outputs.\n",
    "\n",
    "<sup>$^{\\ddagger}$All temporal and spatial inputs/outputs in `Dhara` are stored in `NetCDF4` format.</sup>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1.1 Understanding model inputs\n",
    "##### Import python modules\n",
    "The first thing we'll do is to import all python modules/packages we need for loading and analyzing input data. Modules are generally imported at the very beginning of each notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Numerical packages\n",
    "import numpy as np\n",
    "\n",
    "# NetCDF data format\n",
    "from netCDF4 import Dataset\n",
    "\n",
    "# Plotting package\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LightSource\n",
    "from ipywidgets import *\n",
    "\n",
    "# Seaborn: useful for graphics & statistics\n",
    "import seaborn as sns\n",
    "sns.set_style(style='white')\n",
    "\n",
    "# Bokeh: package for interactive plotting\n",
    "import bokeh.io\n",
    "import bokeh.mpl\n",
    "import bokeh.plotting\n",
    "\n",
    "# functions to make matplotlib and bokeh inline;\n",
    "%matplotlib inline\n",
    "bokeh.io.output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Configuration files\n",
    "We will examine configuration files that are used to navigate the model. Go back to your jupyter main tab, navigate to `$Dhara/examples/1-EBI/` folder and open `ebi.cfg` file.\n",
    "\n",
    "__*TO DO:*__ Open all other configurations file shown in `ebi.cfg` file. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "##### Forcing data\n",
    "Now we know which `netcdf` files will be loaded into Dhara and for what reasons. Specifically, they are `forcings_7536.nc`, `topography.nc`, and 3 vegetation files (`corn.nc`, `miscanthus.nc`, and `switchgrass.nc`).\n",
    "\n",
    "Next, let's examine forcing data which is one of the main drivers of the model. This data is stored in `forcings_7536.nc` as mentioned above.\n",
    "The `plot_1Dinput()` function written below can be used to plot forcing data in jupyter notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Function for loading and plotting 1D forcing input\n",
    "\"\"\"\n",
    "def plot_1Dinput(file_name, var_name):\n",
    "    nc_fid = Dataset(file_name, 'r')\n",
    "    x = nc_fid.variables['DecimalDayOfYear'][:] # common variable on x-axis\n",
    "    var = nc_fid.variables[var_name][:]\n",
    "    \n",
    "    # Plotting data\n",
    "    fig, ax = plt.subplots(figsize=(9,5))\n",
    "    ax.plot(x,var)\n",
    "    ax.set_xlabel('Time (day)')\n",
    "    ax.set_ylabel(var_name)\n",
    "    \n",
    "    # make the plot interactive with Bokeh package\n",
    "    bokeh.plotting.show(bokeh.mpl.to_bokeh())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we want to pass the variable names to an interactive function, so it will allow us to change the variable we want to plot without running the code again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "varname = 'AirTemperature'\n",
    "filename = './examples/1-EBI/forcings_7536.nc' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot data in regular mode\n",
    "plot_1Dinput()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot data in interactive mode\n",
    "# used fixed() if you do not want to change the input\n",
    "interact()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, there are a few more variables in the input forcing file. <br/>\n",
    "__TODO__: Add the following to `varname` and re-plot the forcing data:\n",
    "\n",
    "    * GlobalRadiation\n",
    "    * LongwaveDownward\n",
    "    * Precipitation\n",
    "    * VaporPressureAir\n",
    "    * WindVelocity\n",
    "    * ZenithAngle\n",
    "    * AtmosphericPressure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Vegetation\n",
    "Vegetation netcdf files contain information about Leaf Area Index (`LAI`), Leaf Area Density (`LAD`), and root fraction (`rootfr`). Let's now take a look at these files for the different types of vegetation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Leaf Area Index\n",
    "The first piece of information is LAI. For better comparison, we will plot LAI for corn, miscanthus and switchgrass together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Function for loading and plotting LAI data\n",
    "\"\"\"\n",
    "def plot_LAI(filename):\n",
    "    num_plants = len(filename)\n",
    "    fig, ax = plt.subplots(figsize=(9,5))        \n",
    "    for i in range(0,num_plants):\n",
    "        nc_fid = Dataset(filename[i], 'r')\n",
    "        LAI = nc_fid.variables['LeafAreaIndex'][:]\n",
    "        ax.plot(LAI, label=filename[i])\n",
    "\n",
    "    ax.set_xlabel('Time')\n",
    "    ax.set_ylabel('LAI')\n",
    "    #ax.legend()\n",
    "    bokeh.plotting.show(bokeh.mpl.to_bokeh())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cropnames = ['./examples/1-EBI/corn.nc',\n",
    "             './examples/1-EBI/miscanthus.nc',\n",
    "             './examples/1-EBI/switchgrass.nc']\n",
    "\n",
    "# TODO -- complete the arguments\n",
    "plot_LAI(cropnames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Leaf Area Density and Root fraction\n",
    "\n",
    "Leaf Area Density and root fraction are also important vegetation parameters. This time, we will plot these two parameters side by side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Function for loading and plotting LAI data\n",
    "\"\"\"\n",
    "def plot_LAD_Root(filename):\n",
    "    num_plants = len(filename)\n",
    "    fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(10,5))        \n",
    "    for i in range(0,num_plants):\n",
    "        nc_fid = Dataset(filename[i], 'r')\n",
    "        LAD = nc_fid.variables['LeafAreaDensity'][:]\n",
    "        nl_can = nc_fid.variables['NumCanoyLayers'][:]\n",
    "        ylad = np.linspace(1,nl_can,nl_can)\n",
    "        rootfr = nc_fid.variables['RootFraction'][:]\n",
    "        nl_soil = nc_fid.variables['NumSoilLayers'][:]\n",
    "        ysoil = np.linspace(1,nl_soil,nl_soil)\n",
    "        \n",
    "        ax[0].plot(LAD, ylad/nl_can, label=filename[i])\n",
    "        ax[1].plot(rootfr, -ysoil/nl_soil, label=filename[i])\n",
    "\n",
    "    ax[0].set_xlabel('Leaf Area Density')\n",
    "    ax[0].set_ylabel('z/h [-]')\n",
    "    ax[0].grid(True)\n",
    "\n",
    "    ax[1].set_xlabel('Root fraction')\n",
    "    ax[1].set_ylabel('z/d [-]')\n",
    "    ax[1].grid(True)\n",
    "    ax[1].legend(loc=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# TODO -- complete the arguments\n",
    "plot_LAD_Root(cropnames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Topography\n",
    "We next look at topographic data of the EBI site that will be used for integrated flow model. \n",
    "As we know, topography file of EBI named `topographyEBI.nc` is also in `netCDF` format. \n",
    "The variable inside this netCDF file that stores information on topography is named as `Topography`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_topography(filetopo):\n",
    "    # Load topography in netCDF format\n",
    "    nc_fid = Dataset(filetopo, 'r')\n",
    "    topo = nc_fid.variables['Topography'][:] \n",
    "\n",
    "    # Set font size and light source\n",
    "    font_size = 14\n",
    "    ls = LightSource(azdeg=315, altdeg=45)\n",
    "\n",
    "    # Create a figure with 2 subplots\n",
    "    fig,ax = plt.subplots(ncols=2, nrows=1, figsize=(12,8))\n",
    "\n",
    "\n",
    "    # SUBPLOT 0 ---------\n",
    "    im0 = ax[0].imshow(topo, cmap=plt.cm.Spectral_r) # ... Plot topography\n",
    "    ct = ax[0].contour(topo,\n",
    "                       np.arange(np.floor(topo.min()),np.ceil(topo.max()),1.0),\n",
    "                       linewidths=0.5,\n",
    "                       cmap=plt.cm.jet_r) # ... Plot contour\n",
    "    fig.colorbar(im0, ax=ax[0], pad = 0.1, orientation='horizontal')\n",
    "    ax[0].clabel(ct, inline=True, fmt='%1.1f', fontsize=10) \n",
    "\n",
    "\n",
    "    # SUBPLOT 1 ---------\n",
    "    im1 = ax[1].imshow(ls.shade(topo, cmap=plt.cm.gray_r)) # ... Plot topography \n",
    "                                                           #     with hillshade light source\n",
    "    fig.colorbar(im1,ax=ax[1], pad = 0.1, orientation='horizontal')\n",
    "\n",
    "\n",
    "    # Set axis labels and font size\n",
    "    for i in range(0,2):\n",
    "        ax[i].set_xlabel('X [m]', size=font_size)\n",
    "        ax[i].set_ylabel('Y [m]', size=font_size)\n",
    "        ax[i].tick_params(axis='x', labelsize=font_size)\n",
    "        ax[i].tick_params(axis='y', labelsize=font_size)\n",
    "\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# TODO -- locate topography.nc file and visualize data\n",
    "filetopo = ''\n",
    "plot_topography(filetopo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1.2 Running the code for EBI\n",
    "We have analyzed all configurations and inputs for EBI site. \n",
    "Next, we want to run it on ROGER which has hybrid CPU-GPU parallel computing capability.\n",
    "We will request 16 CPUs & 1 GPU to run simulation as illustrated as below:\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "<img src=\"http://hydrocomplexity.net/images/dhara/ebi-mpi-cuda.png?cache=none\" alt=\"\" width=\"500\">\n",
    "\n",
    "</br>\n",
    "<left>*__Figure 5. MPI-CUDA hybrid computation in the EBI example__ - 16 MPI processes are used for MLCan model on the top and 1 GPU device is used for integrated surface-subsurface flow model. Data are exchanged to couple the two models.*</left>\n",
    "\n",
    "Running computational jobs on the head node of supuercomputers is strictly prohibited. \n",
    "We need to request a compute node or submit a batch job on the queue for computation.\n",
    "The following work must be done on `ROGER` using `terminal`, not `jupterhub`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's run an <span style=\"color:#cc7a00\">__interactive__</span> simulation. After you login to `ROGER`, type the following to request a compute node from the login (or head) node:\n",
    "\n",
    "<pre>\n",
    "    cg-gpu01$ qsub -I -l walltime=00:30:00,nodes=1:ppn=20 -A czomi\n",
    "    cg-gpuXX$ cd $Dhara/examples/1-EBI\n",
    "    cg-gpuXX$ mpirun –np num_proc ../../bin/dhara -option\n",
    "</pre>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we submit a <span style=\"color:#cc7a00\">__batch job__</span> to the system. This work is done by submitting a *portable batch system* file that instructs ROGER to do all the work we want. Open `ebi.pbs`:\n",
    "\n",
    "<pre>\n",
    "    cg-gpu01$ cd $Dhara/examples/1-EBI\n",
    "    cg-gpu01$ nano ebi.pbs\n",
    "</pre>\n",
    "\n",
    "After you input your e-mail address and modify all other information as needed, we submit the job as follow:\n",
    "<pre>\n",
    "    cg-gpu01$ qsub ebi.pbs\n",
    "</pre>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### 3.1.3 Analyzing Model Results\n",
    "\n",
    "Once your simulations are complete, we will use tools similar to those discussed above to analyze the results. Remember that all the results are saved in `netcdf` format as well.\n",
    "\n",
    "##### 2D overland flow visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_overland(time):\n",
    "    nc_f = '../examples/1-EBI/results/olf2D_'+str(time)+'.nc'\n",
    "    nc_fid = Dataset(nc_f, 'r')\n",
    "    water_depth = nc_fid.variables['water_depth'][:]\n",
    "    water_elevation = nc_fid.variables['water_elevation'][:]    \n",
    "    \n",
    "    # Set font size and light source\n",
    "    font_size = 14\n",
    "\n",
    "    # Create a figure with 2 subplots\n",
    "    fig,ax = plt.subplots(ncols=2, nrows=1, figsize=(12,8))\n",
    "\n",
    "    # SUBPLOT 0 ---------\n",
    "    im0 = ax[0].imshow(water_elevation, \n",
    "                       cmap=plt.cm.Spectral_r,\n",
    "                       vmin=water_elevation.min(), \n",
    "                       vmax=water_elevation.max())\n",
    "    fig.colorbar(im0, ax=ax[0], pad = 0.1, orientation='horizontal')\n",
    "    ax[0].set_title('Water Elevation',fontsize=font_size+2)\n",
    "\n",
    "    # SUBPLOT 1 ---------\n",
    "    im1 = ax[1].imshow(water_depth, \n",
    "                       cmap=plt.cm.jet_r,\n",
    "                       vmin=0.0, vmax=0.2)  \n",
    "    cbar = fig.colorbar(im1, ax=ax[1], pad = 0.1, orientation='horizontal')\n",
    "    ax[1].set_title('Water Depth',fontsize=font_size+2)\n",
    "\n",
    "    # Set axis labels and font size\n",
    "    for i in range(0,2):\n",
    "        ax[i].set_xlabel('X [m]', size=font_size)\n",
    "        ax[i].set_ylabel('Y [m]', size=font_size)\n",
    "        ax[i].tick_params(axis='x', labelsize=font_size)\n",
    "        ax[i].tick_params(axis='y', labelsize=font_size)\n",
    "\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# TODO -- visualize overland flow using interact() with: time=(0,7500,100)\n",
    "interact();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### 3D overland flow visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_ssf(time,var,layer):\n",
    "    nc_f = '../examples/1-EBI/results/ssf3D_'+str(time)+'.nc'\n",
    "    nc_fid = Dataset(nc_f, 'r')\n",
    "    data = nc_fid.variables[var][:]\n",
    "    if (var=='moisture'):\n",
    "        vmin=0.3\n",
    "        vmax=0.45\n",
    "        cmap=plt.cm.Spectral\n",
    "    else:\n",
    "        vmin=-1.0\n",
    "        vmax=0.0\n",
    "        cmap=plt.cm.jet_r\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(10,8))\n",
    "    im = ax.imshow(data[layer,:,:], interpolation='bilinear', \n",
    "                                    aspect='auto',\n",
    "                                    vmin=vmin, vmax=vmax,\n",
    "                                    cmap=cmap)\n",
    "    fig.colorbar(im,ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "varssf = ['moisture','pressure_head']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# TODO -- visualize overland flow using interact() with: \n",
    "#    time=(0,7500,100)\n",
    "#    var=varssf\n",
    "#    layer=(0,14,1)\n",
    "interact()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also visualize vertical cross sections of the soil. Make it as an exercise if you want."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Temporal analyses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_2dtime(var):\n",
    "    nc_f = '../examples/1-EBI/results/stat2D.nc'\n",
    "    nc_fid = Dataset(nc_f, 'r')\n",
    "    data = nc_fid.variables[var][:]\n",
    "    if (var=='moisture_mean'):\n",
    "        vmin=0.2\n",
    "        vmax=0.45\n",
    "        cmap=plt.cm.Spectral\n",
    "    else:\n",
    "        vmin=-2.0\n",
    "        vmax=0.0\n",
    "        cmap=plt.cm.jet_r\n",
    "        \n",
    "    fig, ax = plt.subplots(figsize=(18,3))\n",
    "    im = ax.imshow(data.T/len(data), \n",
    "                   interpolation='bilinear', \n",
    "                   vmax=vmax, vmin=vmin,\n",
    "                   aspect='auto',\n",
    "                   cmap=cmap)\n",
    "    fig.colorbar(im,ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "var_data = ['moisture_mean','pressure_mean']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "interact(plot_2dtime, var=var_data);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_1doutput(var_name,proc):\n",
    "    nc_fid = Dataset('../examples/1-EBI/results/proc'+str(proc)+'_mlcan_1d.nc', 'r')\n",
    "    var = nc_fid.variables[var_name][:]\n",
    "    \n",
    "    # Plotting data\n",
    "    fig, ax = plt.subplots(figsize=(10,6))\n",
    "    ax.plot(var)\n",
    "    ax.set_ylabel(var_name)\n",
    "    ax.set_xlabel('Time')    \n",
    "    ax.set_title('MLCan - Proc ' + str(proc))    \n",
    "    \n",
    "    # Want to make it interactive with Bokeh, uncomment line below\n",
    "    bokeh.plotting.show(bokeh.mpl.to_bokeh())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "var1Doutput = ['An_can', 'LE_can', 'H_can', 'TR_can', 'Rnrad_can']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot data in interactive mode\n",
    "interact(plot_1doutput, var_name=var1Doutput, proc=(0,15,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 3.2 Asking questions\n",
    "Now, we have understood how Dhara model work. Let's use it as a tool to ask and explore some interesting scientific questions. For instance:\n",
    "* What are the impacts of elevated CO<sub>2</sub> on ecohydrologic dynamics? \n",
    "* What are the impacts of air temperature increase on evapotranspiration and soil moisture dynamics?\n",
    "* What are their simultaneous impacts on ecohydrology?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Bring your research questions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also important to use Dhara as a framework to study questions related to your research."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
